function [vahat,vapost,ll] = vaposterior(vadist,vadata)

if ~(size(vadata,1)==1)
    error('wrong dimension of value added data')
elseif range([vadata.teacherid])~=0
    error('value added data corresponds to different teachers')
end
    
if isfield(vadist,'residvar')
    Y = vadata.Y;
    Z = vadata.Z;
    ERRCOV = vadist.residvar*eye(size(Y,1));
else
    if size(vadist.residcov,2)~=size(vadata,2)
        error('incorrect multiple outcome data')
    end
    Y = cat(1,vadata.Y);
    Z = blkdiag(vadata.Z);
    uni_sid = unique(cat(1,vadata.sid));
    sid_ind = cellfun(@(sid) ismember(uni_sid,sid),{vadata.sid},'unif',0);
    I = cat(2,sid_ind{:});
    ERRCOV = kron(vadist.residcov,eye(size(I,1)));
    ERRCOV = ERRCOV(I,I);
end

E = vadist.mean;
V = vadist.cov;
pr = vadist.pr;

if (numel(Y)>0)

    ll = zeros(size(pr));
    for c = 1:numel(pr)
        [E(:,c),V(:,:,c),~,ll(c)] = mvnfilter(Y,Z,E(:,c),V(:,:,c),ERRCOV);
    end
    ll = ll + log(pr);
    maxll = max(ll);
    ll = ll - maxll;
    pr = exp(ll)./sum(exp(ll));
    ll = log(sum(exp(ll))) + maxll;

else

    ll = 0;

end

vahat = E*pr';
vapost = struct('mean',E,'cov',V,'pr',pr);
